Orientation control by flow: Exact results and Langevin simulations 
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Rod shaped objects suspended in a flowing liquid might be orientated by the velocity, nature 
of the liquid, the flow and the geometry of the channel containing the flow. Orientation settings 
might enhance or inhibit certain chemical reactions between the objects, other chemicals or with the 
walls of vessels holding the flowing suspension. The probability density function (PDF) describing 
the orientations of rod shaped objects in a flowing liquid satisfies a Fokker-Planck equation whose 
solution is obtained analytically as well as numerically from Langevin simulations for different 
flow parameters. The analytical and numerical methods developed in the present work enable us 
to calculate accurately the PDF for a range of the Peclet number a covering several orders of 
magnitude, 10"'* < a < 10®. We apply these results to the experimental determination of dichroism 
and birefringence of the suspension as a function of a. 
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I. INTRODUCTION 

In many physical, chemical, biological processes, the 
behavior and orientation of rod shaped objects (RSO) 
such as fibers, nanotubes, nanowires, DNA, macro- 
molecules... suspended in a flowing liquid affect the 
transport, rheology, chemical and hydrodynamic char- 
acteristics of the suspension [ij. Orientation control for 
the sake of aligning or separating RSO, enhancing reac- 
tion between them, with other chemicals or vessel walls is 
important in many areas of science and technology such 
as sedimentation, blood flow, pulp and paper, polymer 
processing, microfluidic devices, ferrofluids... 

Objects suspended in flow undergo two types of mo- 
tion: smooth motion due to the average fluid velocity 
field and erratic random motion produced by the fluctu- 
ating fluid velocity, temperature and inertia driven mo- 
tion. The resulting change in the suspension microstruc- 
ture can have a significant effect on the mechanical, ther- 
mal, optical, electrical, magnetic and chemical proper- 
ties. 

The objective of this work is to investigate the effect 
of non-turbulent flow on the rotational diffusion of a di- 
lute suspension of non-interacting RSO in a planar con- 
traction and how its impact might be measured with an 
optical measurement such as dichroism or birefringence 
of the suspension. 

Generally a Fokker-Planck equation accurately models 
the orientation state of non-interacting RSO in hydro- 
dynamic nonhomogcnous flow and the main orientation 
parameter is the rotational Peclet number a that rep- 
resents the interplay between the randomizing effect of 
temperature induced rotational diffusion and the orient- 
ing effect of streamwise mean rate of strain due to the 
flow [1]. 

Two types of flow are of particular importance: simple 
shear flow and elongational (or, extensional) flow. Simple 
shear flow is a velocity profile where the gradient in the 



fluid flow velocity is constant, whereas for elongational 
flow the sample is compressed in one direction and elon- 
gated in the other direction. Such ffows can be either 
stationary or oscillatory. 
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FIG. 1: Geometry of the shear flow and the RSO of length ^ 
making an angle Q with the flow direction along x. The shear 



rate of the flow is 7 = 



The rotational Peclet number measuring the relative 
strength of hydrodynamic interactions and Brownian 
forces is defined mathematically as the ratio of hydro- 
dynamic shear flow and rotational diffusion constant: 



a 



7 



(1) 



the 



where 7 is the hydrodynamic shear rate and Drot 
thermal diffusion 

The flow shear rate defined by 7 = with Vr^ the 
velocity field in the x direction of the flow (see fig. [T]) is 
related to translational degrees of freedom whereas Drot, 
the Brownian diffusion coefficient, governs the rotational 
motion of the RSO around its center of mass and hence 
relates to rotational degrees of freedom. 

The rotational diffusion coefficient, Drot is given by Q|: 



(2) 
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with fcs the Boltzmann constant, T the temperature, 
expresses the influence of flow on orientation. Usually it 
is determined experimentally from RSO orientation stud- 
ies in a constriction with flat walls, t is the RSO length 
d its diameter and -qg the solvent viscosity. 

Non-interacting dilute concentration of RSO of density 
p, yield the conditions: p^^d < 1 and p£^ < 1. 

The Reynolds number defined as [5j]: Re = ^'^"^ is 
much smaller than 2000 in order to have a non-turbulent 
flow, ps is solvent density and L the width of the planar 
contraction (see fig. [T]) . 

Boeder [6| is the first to have studied this problem in 
the bulk of a flowing liquid from the theoretical point of 
view. He derived, without the use of the Langevin for- 
malism or the Fokker-Planck equation Q, an ordinary 
differential equation (ODE) governing the PDF P{9) de- 
scribing the average orientations of the RSO: 



^P{e) + ^[asm\e)Pm=0 (3) 



The above ODE is derived for the motion of RSO, of 
negligible cross-sectional area, in the plane of the flow, 
without any boundary conditions and in the dilute case. 
The angle 9 describes the orientations of the RSO with 
respect to the liquid flow direction (see fig. [T] 

Since the PDF P{0) is 7r-periodic one might write a 
Fourier series solution valid for small and large values 
of a. For small values of a we perform a perturbation 
analysis of the Fourier series whereas for large values of 
a asymptotic analysis might be done [6|] . Many improve- 
ments have been made since, to remove restrictions on 
the cross-sectional areas of the RSO, consider rotational 
diffusion in three dimensions and introduce internal vi- 
brational and rotational degrees of freedom within the 
RSO. 

The purpose of this work is to provide a solution in 
closed form for the above ODE with different analytical 
and numerical methods to obtain P{0) for a wide range of 
the Peclet number a and gauge the impact on orientation 
control and the measurable optical properties. Another 
goal is to compare the analytical approach with Langevin 
simulations of the PDF, under the same conditions in 
the bulk of a flowing liquid. This comparison is useful 
because it can confirm on one hand the robustness of 
the analytical methods. On the other hand it provides a 
necessary limiting bulk condition for similar simulations 
near solid surfaces. 

This paper is organised as follows: In Section II, we 
present an exact analysis of the ODE (see eq. [3]) to ob- 
tain the probability distribution function (PDF), for a 
wide range of a. In Section III, Langevin simulation and 
comparison with the exact results are presented. In sec- 
tion IV, large values of a are treated along with scaling 
results, optical properties and fiow control. The conclu- 
sions are given in Section V. 



II. ANALYTICAL SOLUTION 

In this section, a procedure for the accurate numerical 
analysis of the ODE (see eq. [3]) and its associated prob- 
ability distribution function, P{0), is given for a wide 
range of a. Turbulence effects are known to take place 
for large values of a (typically for a > 10''). Although the 
above ODE (eq.[3]) ceases to apply beyond this limit, the 
present mathematical analysis is not limited by this, and 
numerical solutions may be calculated in our approach 
for values of a ~ 10*. P{0) is the solution of a second- 
order ODE (eq.[3]). The PDF is 7r-periodic since the RSO 
are indistinguishable when oriented at 9 or 9 + tt; hence 
the TT-periodic boundary conditions: 



P(0) = P{tt) and P'(0) = P'(7r) 



(4) 



where P'{9) = ■^P{9). Moreover, the PDF has to be 
normalised over the interval [0, tt]: 



P{9)d9 = 1 



(5) 



The determination of the associated PDF is conse- 
quently a constrained (because of the normalization con- 
dition) boundary value problem (EVP) eq. [3l Neverthe- 
less one may derive two other versions of the ODE eq. [3l 
The first is a ID equation that takes the form: 



P'(e') +asin^(e')P(e') = c 



(6) 



with initial condition: P{9 = 0) = P{0). This is math- 
ematically sound provided the initial value P{0) and the 
constant C are known as shown below. The second form 
is a 2D system that may be presented as: 



^,PiO)= P'iO) 

-t:P'(9) = -asm(9)\sm{9)P'(9) + 2cos(9)P(9)] 
d9 

(7) 

Several technical procedures may be used to determine 
P{9), depending on the different possible formulations of 
the problem, as an initial value problem (eq. [6]and eq.[7]) 
or as a EVP (eq. 2]). In all cases two quantities, the 
constant C in eq. [5] (equal to P'(0) by substituting 9 = 
and assuming finiteness of P{9 — > 0) in eq. O and the 
initial value P(0), have to be evaluated for any value of 
a. We have developed two methods, described below, to 
determine C and P(0). Firstly, a direct method based on 
the solution of eq.[Sl and secondly a minimisation method 
based on a multidimensional secant method (also called 
Broyden method p]). 

The direct method to evaluate C and -P(O) is as follows. 
The formal solution of eq. [Slmay be given generally as: 



3 



P{e) = Cexp[a[sin(26l)/2 - e]/2] x 



(8) 



exp[— a(sin(2a;)/2 — x)/2]dx 



The lower limit — cx) is surprising since the problem is 
defined over the angular interval [0, tt] . It can analytically 
be shown that the lower limit — oo is the only possibility 
compatible with the boundary conditions given by eq. |4l 
Performing a change of variables, the solution may then 
be written as: 



P{e) ^ (2C/a) exp[a sin(26i)/4] x (9) 

/>oo 

/ exp(— x) exp[ a sin(4a;/Q; — 26)/A\dx 
Jo 

The form in eq. [10] is more stable numerically than 
the previous one of eq. [9l since the numerically trouble- 
some exp(0) and exp{—9) terms are avoided. Neverthe- 
less when a increases the exp[a sin(20) /4] term will cause 
problems despite the bounded values of the sine function, 
forcing us to turn to other methods based on extrapola- 
tion techniques. The constant C is next determined from 
the normalisation condition of the PDF (eq. [5]) yielding: 



(f) 



/_\ dej^ dx exp(-x) exp[(f ) sin(2^^) cos(20 - f 



(10) 



whereas P(0) is given by 



P(0) = (2C/a) / exp(-x) exp[ a sin(4a;/a)/4]da; 
Jo 

(11) 

The PDF depends on a which we may want to vary 
over several orders of magnitude. The difficulty in solving 
the ODE stems from the fact that its nature may be 
modified when a increases, turning the problem into a 
singular perturbation one in as detailed in section 
IV. 



III. LANGEVIN SIMULATIONS AND 
COMPARISON WITH THE EXACT EQUATION 

We convert the deterministic equation of motion for 
the RSO into a Langevin equation by adding a random 
force term and analyze statistically the resulting equa- 
tion. 

We start from the angular speed: 



dO 
'dt 



(12) 




7t/2 



FIG. 2: Analytical PDF as a function of 
1. and 10.0 normalised over the interval [- 



? for small a=0.1, 
7r/2,7r/2] 



and add to it a random term X£,{t) with White noise 
statistical properties: 



< ^{t) and < ^{t)^{t') S{t - t') 



(13) 



This means the hydrodynamic forces tend to act on 
the RSO rotating them in the shear flow with an aver- 
age angular speed lo in the presence of perturbating fast 
random orientations represented by £,{t). The equation 
of motion for 9 becomes: 



d9 = -^sin^ {9)dt + XdW{t) 



(14) 



where dW{t) = ^{t)dt is a Wiener process [3]. 

This is an Ito Stochastic Differential Equation (SDE) 
(see for instance Gardiner f^) that can be transformed 
into a ID Fokker-Planck equation that will be identified 
with the ODE equation [3] 

A single stochastic variable 9 ltd SDE {T^ of the form: 



d9 = A{9, t)dt + y^B{9,t)dW{t) (15) 
is equivalent to a ID Fokker-Planck equation: 



- ^[^(^)^(^)] + = (16) 

Identification with the ODE equation[3]yields the value 
of the coefficient A{9) as: A{9) = —^Baam^{9) whereas 
B is a constant independant of 9. Hence the SDE writes: 



d9 



1 



Basm^9)dt + VBdW{t) 



(17) 



The hydrodynamic forces tend to act on the RSO ro- 
tating them in the shear flow with an average angular 
speed w, given by eq. refangspeed. 

This completely characterizes the SDE parameters as: 
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FIG. 3: Langevin simulated PDF as a function of 6 for small 
a=0.1, 1. and 10.0 normalised over the interval [— 7r/2, 7r/2]. 
The analytical curves are displayed for comparison. 



7 = 2^"' 



B = 2Dr 



(18) 



The method of integration we use, is a 4th order 
Stochastic Rungc-Kutta routine: 30-3S-2G developed by 
Helfand ^] where 30 is third order, 3S is three stages 2G 
means we need two Gaussian Random Variables per step. 
Order k means the solution agrees with the average Tay- 
lor expansion of the series solution to order h'' where 
h is the integration step. 

In order to assess the validity of the method, we ran 
several tests and compared the results to known analyt- 
ical solutions. 

The results we find for the Langevin simulated PDF 
are displayed in fig. [3] and fig. [5] for small and large val- 
ues of the Peclet number. The results of the Langevin 
simulation converged after ~ lO"^ steps toward the ana- 
lytically determined PDF. 



IV. LARGE VALUES OF THE PECLET 
NUMBER: SCALING AND OPTICAL 
PROPERTIES 

In order to cope with the wide range over which a 
may vary, we classify the various methods for solving 
the problem according to the value of a. 

For moderate a in the range : [10~'*,10^], we can 
proceed either directly from the analytic solution of the 
first-order ODE or once C and P{0) are known, a simple 
ID Runge-Kutta method is used to solve the first-order 
ODE. 



minimum |P(0) — P(7r) 
P{d)dd ~ 



minimum 



II 



(19) 



In the above minimisation problem, the first condition 
comes from periodicity of the PDF, whereas the second 
expresses the normalisation of the PDF. 

When a becomes large i.e. in the range : [10^,10^], 
we calculate C and P(0) by an extrapolation method 
and solve eq. [6l or the system [71 with singular pertur- 
bation integration methods To illustrate the differ- 
ent numerical methods used in our approach to solve the 
Boeder differential equation, and to calculate the PDF, 
some numerical results are presented, for a wide range of 
values of a. In Figs. 2 and 4, these PDF functions are 
depicted, for moderate (q;=0.1,1,10.) and large values 
(a=100., 1000. and 10,000) of a. The PDF results are 
normalised with respect to unity. The scaling results are 
presented next for large values of a. 



A. Scaling results 

Eq. [To] can be simplified for large a by replacing the 
cosine term in the exponential integrand by 1, since it 
will lead to higher order terms in 1/a and reducing the 
double integration appearing in the C denominator into 
a simpler one, namely: 



C = 



1 

as 



20„ 



, exp(— 2a;3/3)cia:] 



(20) 



where Omax is the value of the angle that corresponds 
to the PDF maximum at which the ODE eq. [3] writes: 



a Sin^ (0 max) P (9 max) = C 

Using the definite integral \u\ : 



(21) 



OO 1 



x''-' exp{-fixP)dx = -irTT(-) (22) 
P P 



we get: 



C 



3a 3 



20 (I) ^r(-) 

^^max V 3 7 V 3 / 



(23) 



Similarly the argument of the exponential integrand in 
eq. [TTJ yields 



In order to find C and P(0) by a minimisation method 
based on a multidimensional secant (Broyden method), 
and a 2D Runge-Kutta integration method to solve the 
two first-order ODE system (Eqs. [7]). We have the fol- 
lowing two conditions: 



F(0) = (2C/a3) / exp(-8xV3Q^)rfx 
Jo 



(2C/3a^)(-l r(-) 



(24) 
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FIG. 4: Analytical PDF as a function of 6 for large q=100.0, 
1000. and 10,000. normalised over the interval [— 7r/2,7r/2] 

The normalisation of the PDF in Eq.l^lwhen a is large, 
yields after using eq. [21] 

6maxP(0max) ~ constant (25) 

which is the area under the PDF curve since the distri- 
bution function becoming sharply peaked around 9max , 
leads to a PDF approximately triangular in shape with 
a height P{9max) and a base equal to 29max- Hence, we 
obtain the following leading behaviours: 

emax^a-^'\ P{e^ax)^a^'^ (26) 
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FIG. 5: Langevin simulated PDF as a function of 6 for large 
a= 100.0, 1000. and 10,000. normalised over the interval 
[— 7r/2, 7r/2]. The analytical curves are displayed for compar- 
ison. 

The exponents controlling the asymptotic behaviour 
above are all consistent with respect to each other, be- 
sides we have checked them numerically up to quite large 
values of a < 10^ . 

Recently, Fry et al. [12] made optical studies of Car- 
bon nanotube suspensions in polyisobutylene and aque- 
ous single-stranded DNA solutions and were able to show 



that shear-induced birefringence and dichroism (SABD) 
An', An" are both proportional to a single quantity 
S —< P2(cos 6) > so-called nematic order parameter that 
is zero for an isotropic distribution of orientations and 1 
for perfect alignment. P-zicosO) is the Legendre poly- 
nomial of order 2 and An' is extracted from the phase 
difference between the transmitted and reference signal 
whereas An" is extracted from the attenuation of light 
intensity. 

They found that SABD scale as a'^'^^^ for values of the 
Peclet number as large as 10^° as they approach the sat- 
uration value of 1. 

Therefore we set out to analyse the scaling of S in the 
Boeder case. Our results shown in fig. [S] indicate that 
in our case, the SABD saturate to 1 and that might be 
attributed to the following argument. 

As the Peclet number increases the orientational PDF 
becomes sharply peaked around Omax that scales as as. 
Since, it is required that the PDF should be normalised 
we ought to have the value of S behaving as: 



S = P{e)P2{cOS0)d0 

- / P{0)P2(cosd)d9 

max 

)P2(C0S 

) (27) 

That implies, given eq. [25l that S behaves as 

P2{cOsd„iax) = ^{3C0S^ 9,nax - 1) 1 - fa^i cloSC tO 

the behaviour found in fig. [5] and obtained from full nu- 
merical integration. 
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FIG. 6: Average value of P2{cos6) versus Peclet number a 
displaying the scaling behaviour at large values of a as it is 
approaching 1. The power law behaviour ~ a^-^^'' for small 
values of a (straight dotted lines) and ~ 1 — 0.93a~'''^* for 
large values of a are indicated. 

Since 9max is a measure of the standard deviation (for 
large values of a) in the fluctuations of the angle 9 about 
the flow main direction, the RSO orientation controlla- 
bility with flow can be analysed with the variation of the 
angle 9max with the Peclet number a. Increasing a makes 



6 



the orientational PDF more sharply peaked around Omax 
and with the requirement that the PDF should be nor- 
malised, Omax approaches the standard deviation. 

In fig. [7] the variation of 9max with the Peclet number 
a is depicted implying that orientation controllability is 
achieved when 9max drops below a small specified value. 
For instance, one has to have a ~ 400 in order to achieve 
Smax < 0.1 radian (see fig. [7]) whereas a '--^ 6 x 10^ is 
required in order to have 9max < 0.01 radian. 



S 0.5 - 
■3 

S 0.4 - 




10 10 ^ lo" 10^ lO'* lo'' 10^ 

a 

FIG. 7: Variation of the angle 6max (in radians) with the 
Peclet number a. For small values of a, 9max saturates at a 
value of about 0.8 radian that is close to -J as seen in fig. [2] 
When a is increased beyond several 100, 9max drops to small 
values indicating that orientation control has set in. 



V. CONCLUSIONS 



The PDF describing the average orientations of RSO in 
a flowing liquid is analytically determined and accurately 
evaluated as a function of a, the rotational Peclet number 
over a wide range. The impact of a on orientation control 
and optical properties is examined. Special analytical as 
well as numerical methods are developed and presented 
in order to calculate accurately this PDF for a wide range 
of a over the interval [10"'* - 10^]. Lan gevin simulations 
agree well with the numerical solutions of the PDF in the 
bulk of the flowing liquid for arbitrary values of a and 
confirm the validity of the analytical solution. Scaling 
results (valid for a > 10^) are also presented. The math- 
ematical nature of the ordinary differential equation that 
the PDF should satisfy is revealed as a singular pertur- 
bation problem when a^^ becomes smaller than about 
10"'^ and special numerical techniques Q ought to be 
used in order to evaluate the PDF and related physical 
quantities. 

The analytical approach, Langevin simulation and scal- 
ing results should be considered equally as complemen- 
tary tools for this type of study and as a valuable guide 
to the harder cases where we have either turbulent flow, 
presence of RSO internal degrees of freedom or more com- 
plex flow geometries. 
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